********************************************************************************
*********************** Do File (6): Tables 5 and 6 ****************************
********************************************************************************

clear
matrix drop _all
set more off

use "$data/EPH_1015_format.dta", clear

*** Keep only observations that have all controls used

keep if age !=. & msa !=. & hhsize !=. & educyr !=. & attsch_ever !=. & native !=. & lit !=. & gender !=. & marstat != . & dec_pcfaminc != .

*** Keep only spouses of domestic workers
keep if sp_domwk != .

********************************************************************************
********************************************************************************

*** Locals for controls

local base_controls "msa year occupation sp_occupation"

local controls "age age2 hhsize lit native attsch_ever educyr educyr2 i.dec_pcfaminc"

*** Label the variables as they appear in the table

label var active "Participation"
label var cont_pens "Registered"
label var hours_mjob "Hours of work per week"
label var inc_mjob_base08 "Income per month"
label var wagehr_mjob_base08 "Wage per hour"

********************************************************************************
********************************************************************************
********************************************************************************

* First run all the regressions to get the p-values and obtain adjusted p-values

qui reghdfe active sp_domwk treat_spdw `controls' if sp_${ctrl_group} == 1, absorb(msa year sp_occupation) vce(cluster msa)
matrix define uPtreat = 2*ttail(e(N),abs(_b[treat_spdw]/_se[treat_spdw]))

qui reghdfe cont_pens sp_domwk treat_spdw `controls' if empstat == 1 & lwagehr_mjob_base08 != . & cont_pens != . & linc_mjob_base08 != . & sp_${ctrl_group} == 1, absorb(`base_controls' `occup') vce(cluster msa)
gen sample_reg = 1 if e(sample)
matrix uPtreat = uPtreat,2*ttail(e(N),abs(_b[treat_spdw]/_se[treat_spdw]))

foreach var in lhours_mjob linc_mjob_base08 lwagehr_mjob_base08 {

	qui reghdfe `var' sp_domwk treat_spdw `controls' if sample_reg == 1, absorb(`base_controls' `occup') vce(cluster msa)
	matrix uPtreat = uPtreat,2*ttail(e(N),abs(_b[treat_spdw]/_se[treat_spdw]))
	
}

* Generate the adjusted p-values and put them in a matrix that will be exported with the regression estimates

matrix uPtreat = uPtreat'
svmat uPtreat, names(unad_p)

qqvalue unad_p1, method(hochberg) qvalue(hochbergP)

qui sum unad_p1
mkmat unad_p1 hochbergP if _n <= r(N), mat(adjPval)
matrix adjPval = adjPval[1..r(N),2]

** Rerun all the regressions, including the adjusted p-values

local i = 1

** Labor force participation

local varlabel: variable label active
qui reghdfe active sp_domwk treat_spdw `controls' if sp_${ctrl_group} == 1, absorb(msa year sp_occupation) vce(cluster msa)
sum active if sp_domwk == 1 & treat == 0 & e(sample)
local meanvar = r(mean)
local qval = round(adjPval[`i',1], 0.001)
outreg2 using "$tables/Table_5a", replace excel keep(treat_spdw) nocons dec(3) label ctitle(`varlabel') addstat(Mean dependent variable,`meanvar', q-value, `qval')
local i = `i' + 1


** Formality conditional on LF participation

local varlabel: variable label cont_pens
qui reghdfe cont_pens sp_domwk treat_spdw `controls' if sample_reg == 1, absorb(`base_controls' `occup') vce(cluster msa)
sum cont_pens if sp_domwk == 1 & treat == 0 & e(sample)
local meanvar = r(mean)
local qval = round(adjPval[`i',1], 0.001)
outreg2 using "$tables/Table_5a", append excel keep(treat_spdw) nocons dec(3) label ctitle(`varlabel') addstat(Mean dependent variable,`meanvar', q-value, `qval')
local i = `i' + 1


foreach var in hours_mjob inc_mjob_base08 wagehr_mjob_base08 {

	local varlabel: variable label `var'
	qui reghdfe l`var' sp_domwk treat_spdw `controls' if sample_reg == 1, absorb(`base_controls' `occup') vce(cluster msa)
	sum `var' if sp_domwk == 1 & treat == 0 & e(sample)
	local meanvar = r(mean)
	local qval = round(adjPval[`i',1], 0.001)
	outreg2 using "$tables/Table_5a", append excel keep(treat_spdw) nocons dec(3) label ctitle(`varlabel') addstat(Mean dependent variable,`meanvar', q-value, `qval')
	local i = `i' + 1
}

********************************************************************************
********************************************************************************

matrix drop _all
use "$data/EPH_1015_format.dta", clear

*** Keep only observations that have all controls used

keep if age !=. & msa !=. & hhsize !=. & native !=. & gender !=. & educyr != .


*** Keep only children aged 16 to 25.

replace par_domwk = . if age > 25
replace par_domwk = . if age < 16

*** Create variable for having completed primary school

gen primary = educlv >= 2
label var primary "Primary school completed"

*** Create variable for having completed secondary school

gen secondary = educlv >= 4
label var secondary "Secondary school completed"

*** Create variable for the years of education of the head of the household
bysort CODUSU nro_hogar time: egen head_educyr = max(cond(relhh == 1, educyr,.))
label var head_educyr "Years of education of the head of the household"

gen head_educyr2 = head_educyr^2

** Treatment for children of domestic workers

gen treat_pardw=treat*par_domwk
label var treat_pardw "Parent is DW x Reform"

keep if par_domwk != .
keep if par_${ctrl_group} == 1

********************************************************************************
********************************************************************************

*** Locals for controls

local base_controls "msa year occupation par_occupation"
local controls "gender age age2 hhsize i.marstat i.dec_pcfaminc head_educyr head_educyr2 singlepar"

*** Label the variables as they appear in the table

label var active "Participation"
label var cont_pens "Registered"
label var hours_mjob "Hours of work per week"
label var inc_mjob_base08 "Income per month"
label var wagehr_mjob_base08 "Wage per hour"

********************************************************************************
********************************************************************************
********************************************************************************



* First run all the regressions to get the p-values and obtain adjusted p-values


*** Labor force participation

qui reghdfe active par_domwk treat_pardw `controls', absorb(msa year par_occupation) vce(cluster msa)
matrix define uPtreat = 2*ttail(e(N),abs(_b[treat_pardw]/_se[treat_pardw]))


*** Formality conditional on LF participation

qui reghdfe cont_pens par_domwk treat_pardw `controls' if lwagehr_mjob_base08 != . & cont_pens != . & linc_mjob_base08 != ., absorb(`base_controls') vce(cluster msa)
gen sample_reg = 1 if e(sample)
matrix uPtreat = uPtreat,2*ttail(e(N),abs(_b[treat_pardw]/_se[treat_pardw]))


foreach var in lhours_mjob linc_mjob_base08 lwagehr_mjob_base08 {

	qui reghdfe `var' par_domwk treat_pardw `controls' if sample_reg == 1, absorb(`base_controls') vce(cluster msa)
	matrix uPtreat = uPtreat,2*ttail(e(N),abs(_b[treat_pardw]/_se[treat_pardw]))

}

********************************************************************************
********************************************************************************

*** Heterogeneity by gender

** Women

*** Labor force participation

qui reghdfe active par_domwk treat_pardw `controls' if gender == 0, absorb(msa year par_occupation) vce(cluster msa)
matrix uPtreat = uPtreat,2*ttail(e(N),abs(_b[treat_pardw]/_se[treat_pardw]))

*** Formality conditional on LF participation

qui reghdfe cont_pens par_domwk treat_pardw `controls' if sample_reg == 1 & gender == 0, absorb(`base_controls') vce(cluster msa)
matrix uPtreat = uPtreat,2*ttail(e(N),abs(_b[treat_pardw]/_se[treat_pardw]))


foreach var in hours_mjob inc_mjob_base08 wagehr_mjob_base08 {

	qui reghdfe l`var' par_domwk treat_pardw `controls' if sample_reg == 1 & gender == 0, absorb(`base_controls') vce(cluster msa)
	matrix uPtreat = uPtreat,2*ttail(e(N),abs(_b[treat_pardw]/_se[treat_pardw]))
	
}


********************************************************************************

** Men

*** Labor force participation

qui reghdfe active par_domwk treat_pardw `controls' if gender == 1, absorb(msa year par_occupation) vce(cluster msa)
matrix uPtreat = uPtreat,2*ttail(e(N),abs(_b[treat_pardw]/_se[treat_pardw]))

*** Formality conditional on LF participation

qui reghdfe cont_pens par_domwk treat_pardw `controls' if sample_reg == 1 & gender == 1, absorb(`base_controls') vce(cluster msa)
matrix uPtreat = uPtreat,2*ttail(e(N),abs(_b[treat_pardw]/_se[treat_pardw]))


foreach var in hours_mjob inc_mjob_base08 wagehr_mjob_base08 {

	qui reghdfe l`var' par_domwk treat_pardw `controls' if sample_reg == 1 & gender == 1, absorb(`base_controls') vce(cluster msa)
	matrix uPtreat = uPtreat,2*ttail(e(N),abs(_b[treat_pardw]/_se[treat_pardw]))

}

********************************************************************************
********************************************************************************

* Generate the adjusted p-values and put them in a matrix that will be exported with the regression estimates

matrix uPtreat = uPtreat'
svmat uPtreat, names(unad_p)

qqvalue unad_p1, method(hochberg) qvalue(hochbergP)

qui sum unad_p1
mkmat unad_p1 hochbergP if _n <= r(N), mat(adjPval)
matrix adjPval = adjPval[1..r(N),2]

** Rerun all the regressions, including the adjusted p-values

local i = 1


*** Labor force participation

local varlabel: variable label active
qui reghdfe active par_domwk treat_pardw `controls', absorb(msa year par_occupation) vce(cluster msa)
sum active if par_domwk == 1 & treat == 0 & e(sample)
local meanvar = r(mean)
local qval = round(adjPval[`i',1], 0.001)
outreg2 using "$tables/Table_5b", replace excel keep(treat_pardw) nocons dec(3) label ctitle(`varlabel') addstat(Mean dependent variable,`meanvar', q-value, `qval') addtext(Controls, Yes, Year Fixed Effects, Yes, Occupation Fixed Effects, No, Metropolitan Area Fixed Effects, Yes, Number of clusters, "`e(N_clust)'")
outreg2 using "$tables/Table_6a", replace excel keep(treat_pardw) nocons dec(3) label ctitle(`varlabel') addstat(Mean dependent variable,`meanvar', q-value, `qval')
local i = `i' + 1


*** Formality conditional on LF participation

local varlabel: variable label cont_pens
qui reghdfe cont_pens par_domwk treat_pardw `controls' if sample_reg == 1, absorb(`base_controls') vce(cluster msa)
sum cont_pens if par_domwk == 1 & treat == 0 & e(sample)
local meanvar = r(mean)
local qval = round(adjPval[`i',1], 0.001)
outreg2 using "$tables/Table_5b", append excel keep(treat_pardw) nocons dec(3) label ctitle(`varlabel') addstat(Mean dependent variable,`meanvar', q-value, `qval') addtext(Controls, Yes, Year Fixed Effects, Yes, Occupation Fixed Effects, Yes, Metropolitan Area Fixed Effects, Yes, Number of clusters, "`e(N_clust)'")
outreg2 using "$tables/Table_6a", append excel keep(treat_pardw) nocons dec(3) label ctitle(`varlabel') addstat(Mean dependent variable,`meanvar', q-value, `qval')
local i = `i' + 1


foreach var in hours_mjob inc_mjob_base08 wagehr_mjob_base08 {

	local varlabel: variable label `var'
	qui reghdfe l`var' par_domwk treat_pardw `controls' if sample_reg == 1, absorb(`base_controls') vce(cluster msa)
	sum `var' if par_domwk == 1 & treat == 0 & e(sample)
	local meanvar = r(mean)
	local qval = round(adjPval[`i',1], 0.001)
	outreg2 using "$tables/Table_5b", append excel keep(treat_pardw) nocons dec(3) label ctitle(`varlabel') addstat(Mean dependent variable,`meanvar', q-value, `qval') addtext(Controls, Yes, Year Fixed Effects, Yes, Occupation Fixed Effects, Yes, Metropolitan Area Fixed Effects, Yes, Number of clusters, "`e(N_clust)'")
	outreg2 using "$tables/Table_6a", append excel keep(treat_pardw) nocons dec(3) label ctitle(`varlabel') addstat(Mean dependent variable,`meanvar', q-value, `qval')
	local i = `i' + 1
}

********************************************************************************
********************************************************************************

*** Heterogeneity by gender

** Women


*** Labor force participation

local varlabel: variable label active
qui reghdfe active par_domwk treat_pardw `controls' if gender == 0, absorb(msa year par_occupation) vce(cluster msa)
sum active if par_domwk == 1 & treat == 0 & e(sample)
local meanvar = r(mean)
local qval = round(adjPval[`i',1], 0.001)
outreg2 using "$tables/Table_6b", replace excel keep(treat_pardw) nocons dec(3) label ctitle(`varlabel') addstat(Mean dependent variable,`meanvar', q-value, `qval')
local i = `i' + 1

*** Formality conditional on LF participation

local varlabel: variable label cont_pens
qui reghdfe cont_pens par_domwk treat_pardw `controls' if sample_reg == 1 & gender == 0, absorb(`base_controls') vce(cluster msa)
sum cont_pens if par_domwk == 1 & treat == 0 & e(sample)
local meanvar = r(mean)
local qval = round(adjPval[`i',1], 0.001)
outreg2 using "$tables/Table_6b", append excel keep(treat_pardw) nocons dec(3) label ctitle(`varlabel') addstat(Mean dependent variable,`meanvar', q-value, `qval')
local i = `i' + 1


foreach var in hours_mjob inc_mjob_base08 wagehr_mjob_base08 {

	local varlabel: variable label `var'
	qui reghdfe l`var' par_domwk treat_pardw `controls' if sample_reg == 1 & gender == 0, absorb(`base_controls') vce(cluster msa)
	sum `var' if par_domwk == 1 & treat == 0 & e(sample)
	local meanvar = r(mean)
	local qval = round(adjPval[`i',1], 0.001)
	outreg2 using "$tables/Table_6b", append excel keep(treat_pardw) nocons dec(3) label ctitle(`varlabel') addstat(Mean dependent variable,`meanvar', q-value, `qval')
	local i = `i' + 1
}


********************************************************************************


** Men

*** Labor force participation

local varlabel: variable label active
qui reghdfe active par_domwk treat_pardw `controls' if gender == 1, absorb(msa year par_occupation) vce(cluster msa)
sum active if par_domwk == 1 & treat == 0 & e(sample)
local meanvar = r(mean)
local qval = round(adjPval[`i',1], 0.001)
outreg2 using "$tables/Table_6c", replace excel keep(treat_pardw) nocons dec(3) label ctitle(`varlabel') addstat(Mean dependent variable,`meanvar', q-value, `qval') addtext(Controls, Yes, Year Fixed Effects, Yes, Occupation Fixed Effects, No, Metropolitan Area Fixed Effects, Yes, Number of clusters, "`e(N_clust)'")
local i = `i' + 1

*** Formality conditional on LF participation

local varlabel: variable label cont_pens
qui reghdfe cont_pens par_domwk treat_pardw `controls' if sample_reg == 1 & gender == 1, absorb(`base_controls') vce(cluster msa)
sum cont_pens if par_domwk == 1 & treat == 0 & e(sample)
local meanvar = r(mean)
local qval = round(adjPval[`i',1], 0.001)
outreg2 using "$tables/Table_6c", append excel keep(treat_pardw) nocons dec(3) label ctitle(`varlabel') addstat(Mean dependent variable,`meanvar', q-value, `qval') addtext(Controls, Yes, Year Fixed Effects, Yes, Occupation Fixed Effects, Yes, Metropolitan Area Fixed Effects, Yes, Number of clusters, "`e(N_clust)'")
local i = `i' + 1


foreach var in hours_mjob inc_mjob_base08 wagehr_mjob_base08 {

	local varlabel: variable label `var'
	qui reghdfe l`var' par_domwk treat_pardw `controls' if sample_reg == 1 & gender == 1, absorb(`base_controls') vce(cluster msa)
	sum `var' if par_domwk == 1 & treat == 0 & e(sample)
	local meanvar = r(mean)
	local qval = round(adjPval[`i',1], 0.001)
	outreg2 using "$tables/Table_6c", append excel keep(treat_pardw) nocons dec(3) label ctitle(`varlabel') addstat(Mean dependent variable,`meanvar', q-value, `qval') addtext(Controls, Yes, Year Fixed Effects, Yes, Occupation Fixed Effects, Yes, Metropolitan Area Fixed Effects, Yes, Number of clusters, "`e(N_clust)'")
	local i = `i' + 1
}